##########################################
########################################## 
## STUDY 2
##########################################
########################################## 
setwd(getwd())
source("libraries.R")

# load data for study 
data_raw <- haven::read_dta(here("Data", "study2.dta"))

##########################################
########################################## 
## Recode Variables
##########################################
########################################## 
df <- data_raw %>% 
  mutate(race = Q35,
         age = Q34,
         edu = Q36,
         female = case_when(
           Q33 == 2 ~ 1,
           Q33 == 1 ~ 0,
           TRUE ~ NA_real_
         ),
         
         #Party Identification (1 = Democrat; 0 = Republican)      
         repdem = case_when(
           Q38_1 == 1 | Q38_1 == 2 | Q38_3 == 1 ~ 0,
           Q38_2 == 1 | Q38_2 == 2 | Q38_3 == 2 ~ 1,
           TRUE ~ NA_real_
         ),
         
         #Dummy indicating whether participant is Independent or not
         independents = case_when(
           Q38_3 == 3 ~ 1,
           TRUE ~ NA_real_
         )) %>% 
  
  #Need for Chaos - Revised
  mutate(need_chaos_r_12 = flip(need_chaos_r_12),
         need_chaos_r_15 = flip(need_chaos_r_15)) %>% 
  mutate(need_chaos = rowMeans(select(., c("need_chaos_1", "need_chaos_2", "need_chaos_3", 
                                           "need_chaos_4", "need_chaos_5", "need_chaos_6", "need_chaos_7"))), 
         na.rm = TRUE) %>% 
  mutate(need_chaos_r = rowMeans(select(., c("need_chaos_r_1", "need_chaos_r_2", "need_chaos_r_3", 
                                             "need_chaos_r_4", "need_chaos_r_5", "need_chaos_r_6", "need_chaos_r_7", "need_chaos_r_12",
                                             "need_chaos_r_15"))), 
         na.rm = TRUE) %>% 
  
  #Political Rumors: Believing and Sharing of DEMocratic and REPublican rumors 
  mutate(true_dem = rowMeans(select(., c("dem1_1","dem2_1","dem3_1"))), na.rm = TRUE) %>% 
  mutate(share_dem = rowMeans(select(., c("dem1_2","dem2_2","dem3_2"))), na.rm = TRUE) %>% 
  mutate(true_rep = rowMeans(select(., c("rep1_1","rep2_1","rep3_1"))), na.rm = TRUE) %>% 
  mutate(share_rep = rowMeans(select(., c("rep1_2","rep2_2","rep3_2"))), na.rm = TRUE) %>% 
  mutate(across(.cols = c(need_chaos, need_chaos_r, female, age, race, edu, true_dem, share_dem, true_rep, share_rep),
                ~ (scale(.) %>%  as.vector),
                .names = "{.col}_zscore")) %>% 
  mutate(across(.cols = c(need_chaos, need_chaos_r, female, age, race, edu, true_dem, share_dem, true_rep, share_rep),
                ~ (zero1(.) %>% as.vector),
                .names = "{.col}01"))

##########################################
########################################## 
## Sample Description 
## (See Supplementary Materials SM1, Table S2)
##########################################
########################################## 
#Gender
print(paste(round(mean(df$female, na.rm=TRUE), 2)*100 , "percent of participants are female"))
#Age
print(paste("Mean age of participants:" , round(mean(df$age, na.rm=TRUE), 2), 
            ", SD:", round(sd(df$age , na.rm =TRUE) , 2)))
#Education
print(paste("Less than high school:", round(prop.table(table(df$edu)), 2)[1]*100, "percent"))
print(paste("High school graduate:", round(prop.table(table(df$edu)), 2)[2]*100, "percent"))
print(paste(" Some college, but no degree:", round(prop.table(table(df$edu)), 2)[3]*100, "percent"))
print(paste("2 year college degree:", round(prop.table(table(df$edu)), 2)[4]*100, "percent"))
print(paste("4 year college degree:", round(prop.table(table(df$edu)), 2)[5]*100, "percent"))
print(paste("Graduate or professional degree:", round(prop.table(table(df$edu)), 2)[6]*100, "percent"))
#Ethnicity
print(paste("White/caucasian:", round(prop.table(table(df$race)), 2)[2]*100, "percent"))

##########################################
########################################## 
## Factor Analysis, etc.: Need for Chaos 
## (See Supplementary Materials SM2)
##########################################
########################################## 
##Exploratory factor analysis: Need for Chaos - Revised (see Table S4)
#All potential items
cfa_efa <- df %>% 
  select(need_chaos_r_1:need_chaos_r_7, need_chaos_r_12:need_chaos_r_17) %>% 
  na.omit()

nfc_revised_efa <- psych::principal(cfa_efa, 
                                    nfactors = 1 , rotate="none" )  
print("Factor loadings:"); nfc_revised_efa$loadings

#Need for Chaos - Revised
cfa_revised <- df %>% 
  select(need_chaos_r_1:need_chaos_r_7 , need_chaos_r_12 , need_chaos_r_15  ) %>% 
  na.omit()
nfc_revised_efa <- psych::principal(cfa_revised, 
                                    nfactors = 1 , rotate="none" )  
print("Factor loadings:"); nfc_revised_efa$loadings

##Need for Chaos - Revised: polychoric correlations
nfc_revised_efa <- psych::principal(cfa_revised, 
                                    nfactors = 1 , rotate="none" , cor = "poly")  
print("Factor loadings:"); nfc_revised_efa$loadings

## Correlation Matrix (See Table S5)
# Note that Item #8 and #11 have been reversed coded compared to Table S5 
cor_matrix <- cor(cfa_revised)
round(cor_matrix, 2)

## Correlation Matrix II (See Table S6)
# Note that Item #8 and #11 have been reversed coded compared to Table S6
polychor(cfa_revised$need_chaos_r_12, cfa_revised$need_chaos_r_1)
polychor(cfa_revised$need_chaos_r_12, cfa_revised$need_chaos_r_2)
polychor(cfa_revised$need_chaos_r_12, cfa_revised$need_chaos_r_3)
polychor(cfa_revised$need_chaos_r_12, cfa_revised$need_chaos_r_4)
polychor(cfa_revised$need_chaos_r_12, cfa_revised$need_chaos_r_5)
polychor(cfa_revised$need_chaos_r_12, cfa_revised$need_chaos_r_6)
polychor(cfa_revised$need_chaos_r_12, cfa_revised$need_chaos_r_7)

polychor(cfa_revised$need_chaos_r_15, cfa_revised$need_chaos_r_1)
polychor(cfa_revised$need_chaos_r_15, cfa_revised$need_chaos_r_2)
polychor(cfa_revised$need_chaos_r_15, cfa_revised$need_chaos_r_3)
polychor(cfa_revised$need_chaos_r_15, cfa_revised$need_chaos_r_4)
polychor(cfa_revised$need_chaos_r_15, cfa_revised$need_chaos_r_5)
polychor(cfa_revised$need_chaos_r_15, cfa_revised$need_chaos_r_6)
polychor(cfa_revised$need_chaos_r_15, cfa_revised$need_chaos_r_7)

#Confirmatory factor analysis: Need for Chaos - Revised (see Table S7)
cfa_revised <- df %>% 
  select(need_chaos_r_1:need_chaos_r_7, need_chaos_r_12, need_chaos_r_15 ) %>% 
  na.omit()

m2  <- ' nfc_oneDim  =~ need_chaos_r_1 + need_chaos_r_2 + need_chaos_r_3 + need_chaos_r_4 + need_chaos_r_5 + need_chaos_r_6 + need_chaos_r_7 + need_chaos_r_12 + need_chaos_r_15  
        need_chaos_r_3~~need_chaos_r_4
        need_chaos_r_3~~need_chaos_r_5
        need_chaos_r_4~~need_chaos_r_5
        need_chaos_r_12~~need_chaos_r_15'

m2 <- lavaan::cfa(m2, data=cfa_revised) 
m2_fit <- summary(m2 , fit.measures=TRUE) 

print(paste("Need for Chaos - Revised,  Chi-square:" , round(m2_fit$FIT["chisq"],2 )))
print(paste("Need for Chaos - Revised, CFI:" , round(m2_fit$FIT["cfi"],2 )))
print(paste("Need for Chaos - Revised, TLI:" , round(m2_fit$FIT["tli"],2 )))
print(paste("Need for Chaos - Revised, RMSEA:" , round(m2_fit$FIT["rmsea"],2 )))

#Cronbach's alpha: Need for Chaos & Need for Chaos - Revised
nfc <- df %>% 
  select(need_chaos_1:need_chaos_7) %>% 
  psych::alpha(.) 
print(paste("Cronbach's alpha, Need for Chaos" , round(nfc$total$raw_alpha,2)))  

nfc_r <- df %>% 
  select(need_chaos_r_1:need_chaos_r_7, need_chaos_r_12, need_chaos_r_15 ) %>% 
  psych::alpha(. , check.keys = TRUE) 
print(paste("Cronbach's alpha, Need for Chaos - revised version:" , round(nfc_r$total$raw_alpha,2)))  

##########################################
########################################## 
## Remove 5% with highest Need for Chaos? (SM8)
## When set to TRUE, used for producing Fig S7, Fig S8, Fig S9
##########################################
##########################################
# Set to TRUE for excluding top 5% of participants with highest Need for Chaos
exclude_obs <- FALSE

if (exclude_obs == TRUE) {
  
  df <- df %>% 
    filter(need_chaos_r_zscore < quantile(need_chaos_r_zscore, probs = .95, na.rm = TRUE))
  
  print("Excluding top 5% observations with highest levels of the Need for Chaos")
  
}

##########################################
########################################## 
## Analysis - Fig. 2 of Main Text
##########################################
########################################## 
#Fig. 2: Right-Hand Panel
hist_df <- df %>% 
  select(need_chaos_r01) %>% 
  na.omit()

mean_nfc_r <- round(mean(hist_df$need_chaos_r01),2)
sd_nfc_r <- sd(hist_df$need_chaos_r01)

fig2_right <- ggplot(hist_df, aes(x=need_chaos_r01)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)), color="grey", fill="grey") + 
  geom_vline(aes(xintercept=mean(need_chaos_r01)),
             color="black",  size=.9) +
  geom_vline(aes(xintercept=mean(need_chaos_r01) + sd(need_chaos_r01)),
             color="black",  size=.9) +
  xlab("Need for Chaos - R") +
  ylab("") +
  theme(axis.text.x = element_text(size=13),
        axis.title = element_text(size=16),
        axis.text.y = element_text(size = 13)) +
  xlim(c(0,1)) + ylim(c(-.02,.15)) +
  annotate("text", x=mean_nfc_r + .075, y=-.015, label="Mean", color = "black", size = 5) +
  annotate("text", x= mean_nfc_r+sd_nfc_r + .076, y=-.015, label="+1 SD", color = "black", size = 5) 

#Save output
save(fig2_right,file=here("Data_recoded","fig2_right.R"))

##########################################
########################################## 
## Analysis - Fig. 3 of Main Text
##########################################
########################################## 
m_share_dem <- lm(share_dem01 ~ need_chaos_r_zscore + repdem + female + edu_zscore + age_zscore + race , data = df )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_r_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study2")

m_share_rep <- lm(share_rep01 ~ need_chaos_r_zscore + repdem + female + edu_zscore + age_zscore + race , data = df )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_r_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study2")

m_true_dem <- lm(true_dem01 ~ need_chaos_r_zscore + repdem + female + edu_zscore + age_zscore + race , data = df )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_r_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study2")

m_true_rep <- lm(true_rep01 ~ need_chaos_r_zscore + repdem + female + edu_zscore + age_zscore + race , data = df )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_r_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study2")

study2_coef <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study2_coef,file=here("Data_recoded","study2_fig3.R"))

#regression table for Supplemental Materials (Table S9)
stargazer(m_share_dem , m_share_rep , m_true_dem , m_true_rep , 
          align=TRUE,
          no.space=TRUE,
          column.labels=c("Share - Dem", "Share - Rep" , "True - Dem" , "Share - Rep"),
          covariate.labels = c("Need for Chaos - R (std.)" , "Party ID (1 = Dem.)" , "Women" , "Education" , "Age  (std.)" , "Ethnicity (1 = non-white)" ,  "Intercept"),
          keep.stat = c("n","ser"),
          dep.var.labels.include = FALSE)

##########################################
########################################## 
## Analysis - Fig. 4 of Main Text
##########################################
########################################## 
#Fig 4: Democratic Voters 
df_dem <- df %>% 
  filter(repdem == 1 )

m_share_dem <- lm(share_dem01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_dem )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study2")

m_share_rep <- lm(share_rep01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_dem )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study2")

m_true_dem <- lm(true_dem01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_dem )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study2")

m_true_rep <- lm(true_rep01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_dem )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study2")

study2_coef_dem <- rbind(share_dem, share_rep, true_dem, true_rep)

#Fig 4: Republican Voters 
df_rep <- df %>% 
  filter(repdem == 0 )

m_share_dem <- lm(share_dem01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_rep )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study2")

m_share_rep <- lm(share_rep01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_rep )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study2")

m_true_dem <- lm(true_dem01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_rep )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study2")

m_true_rep <- lm(true_rep01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_rep )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study2")

study2_coef_rep <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study2_coef_dem,file=here("Data_recoded","study2_fig4_dem.R"))
save(study2_coef_rep,file=here("Data_recoded","study2_fig4_rep.R"))

##########################################
########################################## 
## Analysis - Fig S2 of Supplementary Materials
##########################################
########################################## 
df_independents <- df %>% 
  filter(independents == 1 )

m_share_dem <- lm(share_dem01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_independents )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study2")

m_share_rep <- lm(share_rep01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_independents )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study2")

m_true_dem <- lm(true_dem01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_independents )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study2")

m_true_rep <- lm(true_rep01 ~ need_chaos_r_zscore + female + edu_zscore + age_zscore + race , data = df_independents )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_r_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study2")

study2_coef_independents <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study2_coef_independents,file=here("Data_recoded","study2_figS2.R"))


